# Loughney et al., "Tectonic influence on Cenozoic mammal richness and sedimentation history of the Basin and Range, western North America"
# Science Advances
# R code file for bootstrap analyses, code to plot Figures 2, 3, 4A, and 5A
# final published figures processed in Adobe Illustrator CC

# load required packages and functions file
library(dplyr)
source("bootstrap_functions.r")

# load required species occurrence and locality files for the Basin and Range (br), Northern Basin and Range (nb), Central Basin and Range (cb), and Southern Basin and Range (sb)
brlocs <- read.csv("BR_localities.csv", header = TRUE, stringsAsFactors = FALSE)
brspecies <- read.csv("BR_species_occurrences.csv", header = TRUE, stringsAsFactors = FALSE)
nblocs <- read.csv("NB_localities.csv", header = TRUE, stringsAsFactors = FALSE)
nbspecies <- read.csv("NB_species_occurrences.csv", header = TRUE, stringsAsFactors = FALSE)
cblocs <- read.csv("CB_localities.csv", header = TRUE, stringsAsFactors = FALSE)
cbspecies <- read.csv("CB_species_occurrences.csv", header = TRUE, stringsAsFactors = FALSE)
sblocs <- read.csv("SB_localities.csv", header = TRUE, stringsAsFactors = FALSE)
sbspecies <- read.csv("SB_species_occurrences.csv", header = TRUE, stringsAsFactors = FALSE)

# organize and format data
# create occurrence matrices
brJoin <- dplyr::inner_join(brspecies, brlocs, by = "Locality") # BR
nbJoin <- dplyr::inner_join(nbspecies, nblocs, by = "Locality") # NB
cbJoin <- dplyr::inner_join(cbspecies, cblocs, by = "Locality") # CB
sbJoin <- dplyr::inner_join(sbspecies, sblocs, by = "Locality") # SB

# combine Genus and Species columns
brJoin$SpName <- paste(brJoin$Genus, brJoin$Species, sep = ".") # BR
nbJoin$SpName <- paste(nbJoin$Genus, nbJoin$Species, sep = ".") # NB
cbJoin$SpName <- paste(cbJoin$Genus, cbJoin$Species, sep = ".") # CB
sbJoin$SpName <- paste(sbJoin$Genus, sbJoin$Species, sep = ".") # SB

# create empty columns for bins
brJoin$bin <- rep(0, nrow(brJoin)) # BR
nbJoin$bin <- rep(0, nrow(nbJoin)) # NB
cbJoin$bin <- rep(0, nrow(cbJoin)) # CB
sbJoin$bin <- rep(0, nrow(sbJoin)) # SB

# set parameters
baseAge <- 36		# oldest age of occurrences
topAge <- 0			# youngest age of occurrences
nbins <- baseAge*2	# number of time bins
trials <- 10000		# set bootstrap

set.seed(10000) # set seed for analysis
# run the bootstrap on species richness
brAllBoot <- t(replicate(n = trials, richnessBoot(brJoin, numBins = nbins))) # BR
nbAllBoot <- t(replicate(n = trials, richnessBoot(nbJoin, numBins = nbins))) # NB
cbAllBoot <- t(replicate(n = trials, richnessBoot(cbJoin, numBins = nbins))) # CB
sbAllBoot <- t(replicate(n = trials, richnessBoot(sbJoin, numBins = nbins))) # SB

# calculate median and confidence intervals on species richness
brAllRichness <- ciQuantiles(brAllBoot, 0.95) # BR
nbAllRichness <- ciQuantiles(nbAllBoot, 0.95) # NB
cbAllRichness <- ciQuantiles(cbAllBoot, 0.95) # CB
sbAllRichness <- ciQuantiles(sbAllBoot, 0.95) # SB

# run bootstrap, calculate median and confidence intervals on localities
# bootstrap
brLocBoot <- t(replicate(n = trials, localityBoot(brJoin, numBins = nbins))) # BR
nbLocBoot <- t(replicate(n = trials, localityBoot(nbJoin, numBins = nbins))) # NB
cbLocBoot <- t(replicate(n = trials, localityBoot(cbJoin, numBins = nbins))) # CB
sbLocBoot <- t(replicate(n = trials, localityBoot(sbJoin, numBins = nbins))) # SB

# calculate median and confidence intervals on localities
brBinLocs <- ciQuantiles(brLocBoot, 0.95) # BR
nbBinLocs <- ciQuantiles(nbLocBoot, 0.95) # NB
cbBinLocs <- ciQuantiles(cbLocBoot, 0.95) # CB
sbBinLocs <- ciQuantiles(sbLocBoot, 0.95) # SB

# ------------------------------------------------------------------------------------------------------------------------------------
# plot figures from bootstrapped species richness and number of localities

# Figure 2 - deformation rates, area-change rates, sediment-accumulation rates of Macrostrat packages, and species richness of the Basin and Range
# load additional files
macrostrat <- read.csv("macrostrat_June-2021.csv", header = TRUE, stringsAsFactors = FALSE)
area <- read.csv("area_change.csv", header = TRUE, stringsAsFactors = FALSE)
deformation <- read.csv("deformation_rates_January-2021.csv", header = TRUE, stringsAsFactors = FALSE)

# plot
plot.new();
par(mfrow=c(1,1))
par(oma = c(1, 1, 1, 1));
par(mar = c(5, 4, 4, 2));
plot.window(xlim =c(37, -0.5), ylim = c(-1500, 12000));
axis4lab <- seq(0, 240, 40)
axis(side = 1, at = seq(0, 36, by = 2), cex.axis = 0.9, pos = -1000)
mtext("Age (Ma)", side = 1, line = 2)
axis(side = 2, at = seq(0, 12000, by = 2000), cex.axis = 0.9, pos = 37, las = 1)
mtext("Sediment-accumulation rate (m/Myr)", side = 2, line = 3) # SAR & area-change rate
axis(side = 4, at = seq(0, 12000, by = 2000), lab = axis4lab, cex.axis = 0.9, pos = -1, las = 1)
mtext("Deformation rate (km/Myr)", side = 4, line = 3) # deformation rates & mammal richness
# time scale
rect(xleft = c(36, 36), ybottom = c(-1000, -1000), xright = c(33.9, 33.9), ytop = c(-300, -300)) # Eocene
text(x = baseAge-(baseAge-33.9)/2, y = -600, "Eo", cex = 0.9)
rect(xleft = c(33.9, 33.9), ybottom = c(-1000, -1000), xright = c(23.03, 23.03), ytop = c(-300, -300)) # Oligocene
text(x = 33.9-(33.9-23.03)/2, y = -600, "Oligocene", cex = 0.9)
rect(xleft = c(23.03, 23.03), ybottom = c(-1000, -1000), xright = c(5.33, 5.33), ytop = c(-300, -300)) # Miocene
text(x = 23.03-(23.03-5.33)/2, y = -600, "Miocene", cex = 0.9)
rect(xleft = c(5.33, 5.33), ybottom = c(-1000, -1000), xright = c(1.8, 1.8), ytop = c(-300, -300)) # Pliocene
text(x = 5.33-(5.33-1.8)/2, y = -600, "Plio", cex = 0.9)
rect(xleft = c(1.8, 1.8), ybottom = c(-1000, -1000), xright = c(0, 0), ytop = c(-300, -300)) # Pleistocene
text(x = 1.8-(1.8-topAge)/2, y = -600, "Q", cex = 0.9)
#
lines(macrostrat$mid_bin[2:72], macrostrat$BR_SAR[2:72], col = "#66A61E", lwd = 2) # macrostrat units
lines(macrostrat$mid_bin, area$BR_rate_change_sum, col = "#d95f02", lwd = 2)
lines(macrostrat$mid_bin, brall$BR_rate*50, col = "#7570b3", lwd = 2)
lines(macrostrat$mid_bin, brAllRichness[ , 2]*50, col = "black", lwd = 2)
legend(35, 12000, c("Species richness", "SAR", "Deformation", "Area change"), c("black", "#66A61E", "#7570b3", "#d95f02"))

# -----------------------------------------------------------------------------------------------------------------------------
# Figure 3 - number of fossil localities and fossiliferous units
# load additional files
fossils <- read.csv("fossiliferous_units_January-2021.csv", header = TRUE, stringsAsFactors = FALSE)

# plot
plot.new();
par(oma = c(1, 1, 1, 1));
par(mar = c(2, 4, 1, 2));
par(mfrow = c(3, 1));
CBaxis4lab <- c(0, 20, 40, 60, 80, 100)
axis2lab <- c(0, 5, 10, 15, 20, 25)
# NB 
nbAvgRange <- 1.4
plot(fossils$mid_bin, fossils$NB_number, type = "n", xlim = c(36, 0), ylim = c(0, 50), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 50, by = 10), lab = axis2lab, cex.axis = 0.9, pos = 37, las = 1) # No units
mtext("No. Units", side = 2, line = 3, cex.axis = 1)
axis(side = 4, at = seq(0, 50, by = 10), cex.axis = 0.9, pos = -1, las = 1) # No locs
mtext("No. Localities", side = 4, line = 3, cex.axis = 1)
lines(fossils$mid_bin, fossils$NB_number*2, col = "#65A9DD", lwd = 2, lty = 1, lend = "butt")
points(fossils$mid_bin[2:72], nbBinLocs[2:72, 2], pch = 23, bg = "steelblue", cex = 1.8)
segments(x0 = fossils$mid_bin[2:72]+nbAvgRange, y0 = nbBinLocs[2:72, 2], x1 = fossils$mid_bin[2:72]-nbAvgRange, y1 = nbBinLocs[2:72, 2], col = "gray") # age uncertainty
segments(x0 = fossils$mid_bin[2:72], y0 = nbBinLocs[2:72, 1], x1 = fossils$mid_bin[2:72], y1 = nbBinLocs[2:72, 3], col = "gray") # locality CI
text(33, 45, "Northern BR", cex = 1.5)
# CB
cbAvgRange <- 1.7
plot(fossils$mid_bin, fossils$CB_number, type = "n", xlim = c(36, 0), ylim = c(0, 50), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 50, by = 10), lab = axis2lab, cex.axis = 0.9, pos = 37, las = 1) # No units
mtext("No. Units", side = 2, line = 3, cex.axis = 1)
axis(side = 4, at = seq(0, 50, by = 10), lab = CBaxis4lab, cex.axis = 0.9, pos = -1, las = 1) # No locs
mtext("No. Localities", side = 4, line = 3, cex.axis = 1)
lines(fossils$mid_bin, fossils$CB_number*2, col = "#CD4F39", lwd = 2, lty = 1, lend = "butt")
points(fossils$mid_bin[2:72], cbBinLocs[2:72, 2], pch = 23, bg = "tomato", cex = 1.8)
segments(x0 = fossils$mid_bin[2:72]+cbAvgRange, y0 = cbBinLocs[2:72, 2], x1 = fossils$mid_bin[2:72]-cbAvgRange, y1 = cbBinLocs[2:72, 2], col = "gray") # age uncertainty
segments(x0 = fossils$mid_bin[2:72], y0 = cbBinLocs[2:72, 1], x1 = fossils$mid_bin[2:72], y1 = cbBinLocs[2:72, 3], col = "gray") # locality CI
text(33, 45, "Central BR", cex = 1.5)
# SB
sbAvgRange <- 1.7
plot(fossils$mid_bin, fossils$SB_number, type = "n", xlim = c(36, 0), xlab = "Age (Ma)", ylim = c(0, 50), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
mtext("Age (Ma)", side = 1, line = 3, cex.lab = 1)
axis(side = 2, at = seq(0, 50, by = 10), lab = axis2lab, cex.axis = 0.9, pos = 37, las = 1) # No units
mtext("No. Units", side = 2, line = 3, cex.axis = 1)
axis(side = 4, at = seq(0, 50, by = 10), cex.axis = 0.9, pos = -1, las = 1) # No locs
mtext("No. Localities", side = 4, line = 3, cex.axis = 1)
lines(fossils$mid_bin, fossils$SB_number*2, col = "#F0A63E", lwd = 2, lty = 1, lend = "butt")
points(fossils$mid_bin[2:72], sbBinLocs[2:72, 2], pch = 23, bg = "gold2", cex = 1.8)
segments(x0 = fossils$mid_bin[2:72]+sbAvgRange, y0 = sbBinLocs[2:72, 2], x1 = fossils$mid_bin[2:72]-sbAvgRange, y1 = sbBinLocs[2:72, 2], col = "gray") # age uncertainty
segments(x0 = fossils$mid_bin[2:72], y0 = sbBinLocs[2:72, 1], x1 = fossils$mid_bin[2:72], y1 = sbBinLocs[2:72, 3], col = "gray") # locality CI
text(33, 45, "Southern BR", cex = 1.5)

# ----------------------------------------------------------------------------------------------------------------------------------------
# Figure 4 A - mammalian species richness of BR subregions 
plot.new();
par(mfrow = c(1,1))
par(oma = c(1, 1, 1, 1));
par(mar = c(5, 4, 4, 2));
plot.window(xlim = c(37, -0.5), ylim = c(-30, 150));
axis(side = 1, at = seq(0, 36, by = 2), cex.axis = 0.9, pos = -16)
mtext("Age (Ma)", side = 1, line = 1, cex.lab = 1)
axis(side = 2, at = seq(0, 150, by = 25), cex.axis = 0.9, pos = 37, las = 1)
mtext("No. Species", side = 2, line = 2, cex.axis = 1)
# time scale
rect(xleft = c(baseAge, baseAge), ybottom = c(-15, -15), xright = c(33.9, 33.9), ytop = c(-5, -5)) # Eocene
text(x = baseAge-(baseAge-33.9)/2, y = -10, "Eo", cex = 0.9)
rect(xleft = c(33.9, 33.9), ybottom = c(-15, -15), xright = c(23.03, 23.03), ytop = c(-5, -5)) # Oligocene
text(x = 33.9-(33.9-23.03)/2, y = -10, "Oligocene", cex = 0.9)
rect(xleft = c(23.03, 23.03), ybottom = c(-15, -15), xright = c(5.33, 5.33), ytop = c(-5, -5)) # Miocene
text(x = 23.03-(23.03-5.33)/2, y = -10, "Miocene", cex = 0.9)
rect(xleft = c(5.33, 5.33), ybottom = c(-15, -15), xright = c(1.8, 1.8), ytop = c(-5, -5)) # Pliocene
text(x = 5.33-(5.33-1.8)/2, y = -10, "Plio", cex = 0.9)
rect(xleft = c(1.8, 1.8), ybottom = c(-15, -15), xright = c(0, 0), ytop = c(-5, -5)) # Pleistocene
text(x = 1.8-(1.8-0)/2, y = -10, "Q", cex = 0.9)
# richness curves
lines(nbAllRichness[, 4], nbAllRichness[, 2], col = "#65A9DD", lwd = 2, lend = "butt")
polygon(x = c(nbAllRichness[, 4], rev(nbAllRichness[, 4])), y = c(nbAllRichness[, 3], rev(nbAllRichness[, 1])), col = adjustcolor("#65A9DD", alpha.f = 0.30), border = NA)
lines(cbAllRichness[, 4], cbAllRichness[, 2], col = "#CD4F39", lwd = 2, lend = "butt")
polygon(x = c(cbAllRichness[, 4], rev(cbAllRichness[, 4])), y = c(cbAllRichness[, 3], rev(cbAllRichness[, 1])), col = adjustcolor("#CD4F39", alpha.f = 0.30), border = NA)
lines(sbAllRichness[, 4], sbAllRichness[, 2], col = "#FEC124", lwd = 2, lend = "butt")
polygon(x = c(sbAllRichness[, 4], rev(sbAllRichness[, 4])), y = c(sbAllRichness[, 3], rev(sbAllRichness[, 1])), col = adjustcolor("#FEC124", alpha.f = 0.30), border = NA)
legend(34, 150, c("Northern BR", "Central BR", "Southern BR"), c("#65A9DD", "#CD4F39", "#FEC124"))

# ----------------------------------------------------------------------------------------------------------------------------------------
# Figure 5 A - BR Subregions mammalian species richness
plot.new();
par(oma = c(1, 1, 1, 1));
par(mar = c(2, 4, 1, 2));
par(mfrow = c(3, 1));
# NB 
plot(fossils$mid_bin[2:72], nbAllRichness[2:72, 2], type = "n", xlim = c(36, 0), ylim = c(0, 70), las = 1, ylab = "No. species", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 70, by = 10), cex.axis = 0.9, las = 1)
lines(fossils$mid_bin[2:72], nbAllRichness[2:72, 2], col = "#65A9DD", lwd = 2, lty = 1, lend = "butt")
polygon(x = c(fossils$mid_bin[2:72], rev(fossils$mid_bin[2:72])), y = c(nbAllRichness[2:72, 3], rev(nbAllRichness[2:72, 1])), col = adjustcolor("#65A9DD", alpha.f = 0.30), border = NA)
text(33, 65, "Northern BR", cex = 1.5)
# CB
plot(fossils$mid_bin[2:72], cbAllRichness[2:72, 2], type = "n", xlim = c(36, 0), ylim = c(0, 70), las = 1, ylab = "No. species", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 70, by = 10), cex.axis = 0.9, las = 1)
lines(fossils$mid_bin[2:72], cbAllRichness[2:72, 2], col = "#CD4F39", lwd = 2, lty = 1, lend = "butt")
polygon(x = c(fossils$mid_bin[2:72], rev(fossils$mid_bin[2:72])), y = c(cbAllRichness[2:72, 3], rev(cbAllRichness[2:72, 1])), col = adjustcolor("#CD4F39", alpha.f = 0.30), border = NA)
text(33, 65, "Central BR", cex = 1.5)
# SB
plot(fossils$mid_bin[2:72], sbAllRichness[2:72, 2], type = "n", xlim = c(36, 0), xlab = "Age (Ma)", ylim = c(0, 70), las = 1, ylab = "No. species", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 70, by = 10), cex.axis = 0.9, las = 1)
lines(fossils$mid_bin[2:72], sbAllRichness[2:72, 2], col = "#FEC124", lwd = 2, lty = 1, lend = "butt")
polygon(x = c(fossils$mid_bin[2:72], rev(fossils$mid_bin[2:72])), y = c(sbAllRichness[2:72, 3], rev(sbAllRichness[2:72, 1])), col = adjustcolor("#FEC124", alpha.f = 0.30), border = NA)
text(33, 65, "Southern BR", cex = 1.5)